﻿#region Header

#region License

/*
Copyright (c) 2009, G.W. van der Vegt
All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided
that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this list of conditions and the
  following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and
  the following disclaimer in the documentation and/or other materials provided with the distribution.

* Neither the name of G.W. van der Vegt nor the names of its contributors may be
  used to endorse or promote products derived from this software without specific prior written
  permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#endregion License

#region Changelog

//----------   ---   -------------------------------------------------------------------------------
//Purpose:           Sparse Matrix Class (Compressed Row).
//Author:            Ir. G.W. van der Vegt
//Url:               http://swiss.codeplex.com
//License:           New BSD License
//----------   ---   -------------------------------------------------------------------------------
//dd-mm-yyyy - who - description
//----------   ---   -------------------------------------------------------------------------------
//16-01-2010 - veg - Created.
//                 x In contrast to Delphi version, Arrays are zero based!
//                 - Added Clear();
//                 x Added Boolean result to First() to signal Empty Arrays (After a Clear() call);
//                 - Added this[row,col] as substitute for Element[].
//                 - Made Sparse Matrix Class use a Generic <T> as Value type.
//                   Currently only Int32 and Double are supported!
//                 - Added Determinant().
//                 - Got the multiplication of Generics working by using Convert.ChangeType to convert
//                   back and forth between <T> and Double or Int32.
//                 - Added Invert().
//19-01-2010 - veg - Added copy constructor.
//                 - Moved initialize code from Constructor to Init() method.
//                 - Implemented Compress (Removes zero/default values).
//                 - Swapped Convert.ChangeType() for T.Equals().
//                 - Updated Test Methods.
//22-01-2010 - veg - Added Swiss.cs, based on an idea of Roger Alsing (http://www.rogeralsing.com) in 2008.
//                 - Added SparseElement structure.
//                 - Added IEnumerator interface so one can use foreach (Sparse<Double>.SparseElement se in s) {}.
//                 - Fixed bug in First(), screwing up arrays starting with empty rows!
//                 x Added Transpose().
//                 - Added SaveHarwellBoeing().
//27-01-2010 - veg - Added Cofactor().
//                 - Fixed 'bug' in IEnumerator. When getEnumerator is called more than once, Reset() is not called and the 
//                   subsequent foreach loops do not do their work. Workaround is to call Reset() manually in GetEnumerator()
//                   or implement the Enumerator in a separate class.
//                   see: http://bytes.com/topic/c-sharp/answers/245096-foreach-ienumerator
//                 - Updated code of Determinant() to using a foreach loop.
//                 - Made One a public property in the Swiss class and adjusted the Get property implementation accordingly.
//29-01-2010 - veg * Added Inplace Compress().
//                 - Added Dot().
//                 - Added Extract().
//                 - Added Distance().
//                 - Added Cosine().
//                 x Added Frobenius().
//                 - Added Multiply().
//                 x Added Reduce() (Inplace Exctract)
//                 - Added SaveExcel().
//                 - Corrected Escaping of CSV output.
//30-01-2010 - veg * Added LoadHardwelBoeing().
//                   This implementation does not Transpose inside the loading.
//                   Thus the Matrix loaded in memory has exactly the same layout as the file. 
//09-02-2010 - veg * Changed implementation of Get/SetAsArray functions so they use the dimension 
//                   of the Arrays instead of the rows/columns parameters of the Delphi Functions.
//                 - Implemented a first version of SaveNonSparse().
//                 - Added Auto NonSparse Saving to SaveHarwellBoeing().
//                 - Swapped StringBuilder for FileStream/StreamWriter() so memory load is lower.
//                 - Added Some Profiling on Saving Matrices.
//----------   ---   -------------------------------------------------------------------------------
//dd-mm-yyyy         todo's
//                 * Add Sanity Checks (-1 for RowPtr etc).
//                 * Test Cosine().
//                 * Test Fronebius().
//                 * Turn Console Writes into Asserts().
//                 - Debug SaveMatrixMarket().
//----------   ---   -------------------------------------------------------------------------------
//dd-mm-yyyy         notes
//                 - Changes with a 'x' differ from the Borland Delphi version.
//----------   ---   -------------------------------------------------------------------------------

#endregion Changelog

#endregion Header

namespace Swiss
{
    using System;
    using System.Collections.Generic;
    using System.Diagnostics;
    using System.Linq;
    using System.Text;
    using Swiss;
    using System.Collections;
    using System.IO;
    using System.Reflection;
    using System.Collections.Specialized;

    /// <summary>
    /// Creates a Sparse Matrix for the specified DataType T
    /// </summary>
    /// <typeparam name="T">The type of the Values stored in the Sparse Matrix</typeparam>
    public class Sparse<T> : IEnumerator, IDisposable
    {
        #region Fields

        /// <summary>
        /// A Key for Identifying Harwell Boeing Sparse Matrices.
        /// </summary>
        public static String HBkey = "NetLsa";

        #endregion Fields

        #region Constants

        /// <summary>
        /// Number of Values per line of data for Harwell Boeing File. 8 is equal to SVDLIBC.
        /// </summary>
        internal readonly Int32 HBSplit = 8;

        #endregion Constants

        #region Constructors

        /// <summary>
        /// Creates a Zero Element/Zero Dimensions Matrix
        /// </summary>
        public Sparse()
            : this(0, 0, 0, "")
        {
            //Nothing
        }

        /// <summary>
        /// Initializes a Sparse Matrix.
        /// </summary>
        /// <param name="RowCount">Number of Rows</param>
        /// <param name="ColCount">Number of Columns</param>
        /// <param name="NonZeroCount">Number of NonZero Elements</param>
        public Sparse(Int32 RowCount, Int32 ColCount, Int32 NonZeroCount)
            : this(RowCount, ColCount, NonZeroCount, "")
        {
            //Nothing
        }

        /// <summary>
        /// Initializes a Sparse Matrix.
        /// </summary>
        /// <param name="RowCount">Number of Rows</param>
        /// <param name="ColCount">Number of Columns</param>
        /// <param name="NonZeroCount">Number of NonZero Elements</param>
        /// <param name="Comment">Optional Comment to be written in various output format</param>
        public Sparse(Int32 RowCount, Int32 ColCount, Int32 NonZeroCount, String Comment)
        {
            ExcelSeparator = ';';

            Init(RowCount, ColCount, NonZeroCount, Comment);
        }

        /// <summary>
        /// Copy Constructor.
        /// </summary>
        /// <param name="previous"></param>
        public Sparse(Sparse<T> previous)
        {
            this.ColCount = previous.ColCount;
            this.RowCount = previous.RowCount;
            this.NonZeroCount = previous.NonZeroCount;
            this.Comment = previous.Comment;

            this.Values = previous.Values;
            this.ColNdx = previous.ColNdx;
            this.RowPtr = previous.RowPtr;

            FNdx = previous.FNdx;

            FFN = previous.FFN;
            FFR = previous.FFR;
        }
        #endregion Constructors

        #region Properties

        /// <summary>
        /// Returns the Value passed to Init.
        /// </summary>
        public Int32 ColCount
        {
            get;
            set;
        }

        /// <summary>
        /// Literature - JA  Column Index for each Value.
        /// </summary>
        public Int32[] ColNdx
        {
            get { return FColNdx; }
            set { FColNdx = value; }
        }
        private Int32[] FColNdx;

        /// <summary>
        /// Comment to be saved in various matrix output formats.
        /// </summary>
        public String Comment
        {
            get;
            set;
        }

        /// <summary>
        /// Returns the Value passed to Init.
        /// </summary>
        public Double Fill
        {
            get
            {
                if (RowCount * ColCount != 0)
                {
                    return (100 * (NonZeroCount / (RowCount * ColCount)));
                }
                else
                {
                    return 100;
                }
            }
        }

        /// <summary>
        /// Returns the Value passed to Init.
        /// </summary>
        public Int32 NonZeroCount
        {
            get;
            set;
        }

        /// <summary>
        /// Used when walking the array with First/Next Iterator.
        /// </summary>
        public Int32 Position
        {
            get { return FFN; }
            //set { FFN = value; }
        }

        /// <summary>
        /// Returns the Value passed to Init.
        /// </summary>
        public Int32 RowCount
        {
            get;
            set;
        }

        /// <summary>
        /// Literature - IA - Points to start of Rows in ColNdx.
        /// </summary>
        public Int32[] RowPtr
        {
            get { return FRowPtr; }
            set { FRowPtr = value; }
        }
        private Int32[] FRowPtr;

        /// <summary>
        /// Percentage of Zero Values.
        /// </summary>
        public Double SparseNess
        {
            get { return 100 - Fill; }
        }

        /// <summary>
        ///Literature - A - The actual Sparse Data.
        /// </summary>
        public T[] Values
        {
            get { return FValues; }
            set { FValues = value; }
        }
        private T[] FValues;

        /// <summary>
        /// Used when walking the array with First()/Next() Iterator.
        /// </summary>
        public Int32 FFR
        {
            get;
            set;
        }

        /// <summary>
        /// Used when walking the array with First()/Next() Iterator.
        /// </summary>
        public Int32 FFN
        {
            get;
            set;
        }

        /// <summary>
        /// Used when filling the Matrix with AddElement().
        /// </summary>
        internal Int32 FNdx
        {
            get;
            set;
        }

        /// <summary>
        /// Used in SaveExcel.
        /// </summary>
        internal Char ExcelSeparator
        {
            get;
            set;
        }

        /// <summary>
        /// String.Format code used to format values in textual output files.
        /// </summary>
        public String SparseFmt
        {
            get
            {
                if (String.IsNullOrEmpty(FSparseFmt))
                {
                    SparseFmt = "{0:0.0000}";
                }
                return FSparseFmt;
            }
            set
            {
                FSparseFmt = value;
            }
        }
        private String FSparseFmt = null;

        #endregion Properties

        #region Indexers

        /// <summary>
        /// Allows direct addressing of Elements.
        /// 
        /// Note: Writing can only be done at NonZero Elements!
        /// </summary>
        /// <param name="Row">The Row Index</param>
        /// <param name="Col">The Columns Index</param>
        /// <returns>The Value</returns>
        public T this[Int32 Row, Int32 Col]
        {
            get
            {
                T Result = default(T);

                Int32 e = -1;
                if (Row < RowPtr.Length)
                {
                    e = RowPtr[Row + 1] - RowPtr[Row];
                }

                for (Int32 i = 0; i < e; i++)
                {
                    if (ColNdx[RowPtr[Row] + i] == Col)
                    {
                        Result = Values[RowPtr[Row] + i];

                        break;
                    }
                }

                return Result;
            }

            set
            {
                Int32 e = -1;
                if (Row < RowPtr.Length)
                {
                    e = RowPtr[Row + 1] - RowPtr[Row];
                }

                for (Int32 i = 0; i < e; i++)
                {
                    if (ColNdx[RowPtr[Row] + i] == Col)
                    {
                        Values[RowPtr[Row] + i] = value;

                        break;
                    }
                }
            }
        }

        #endregion Indexers

        #region Methods

        /// <summary>
        /// Adds an Element to the Sparse Matrix.
        /// 
        /// Note: Elements must be added in the correct order (Row by Row) 
        ///       and the number of calls should be equal to the NonZeroCount.
        /// </summary>
        /// <param name="Row">The Row Index</param>
        /// <param name="Column">The Column Index</param>
        /// <param name="Element">The Value</param>
        public void AddElement(Int32 Row, Int32 Column, T Element)
        {
            Values[FNdx] = Element;

            ColNdx[FNdx] = Column;

            if (RowPtr[Row] == -1)
            {
                RowPtr[Row] = FNdx;
            }
            FNdx++;

            if (FNdx == Values.Length)
            {
                //Fill Empty Rows
                for (Int32 i = RowPtr.Length - 2; i >= 0; --i)
                {
                    //if (RowPtr[i] == 0)
                    //{
                    //    break;
                    //}
                    if (RowPtr[i] == -1)
                    {
                        RowPtr[i] = RowPtr[i + 1];
                    }
                }
            }
        }

        /// <summary>
        /// Adds an Element to the Sparse Matrix.
        /// 
        /// Note: Elements must be added in the correct order (Row by Row) 
        ///       and the number of calls should be equal to the NonZeroCount.
        /// </summary>
        /// <param name="se">The Element to Add</param>
        public void AddElement(SparseElement se)
        {
            AddElement(se.Row, se.Column, se.Value);
        }

        /// <summary>
        /// Clears the Matrix and sets NonZeroElements to 0.
        /// </summary>
        public void Clear()
        {
            this.ColCount = 0;
            this.RowCount = 0;
            this.NonZeroCount = 0;

            this.Comment = "";

            this.Values = new T[NonZeroCount];
            this.ColNdx = new Int32[NonZeroCount];
            this.RowPtr = new Int32[RowCount + 1];

            RowPtr[RowCount] = NonZeroCount;

            FNdx = 0;

            FFN = 0;
            FFR = 0;
        }

        /// <summary>
        /// Calculates the Determinant.
        /// 
        /// Note: The Method used ONLY supports Diagonal Matrices.
        /// </summary>
        /// <returns>The Determinant</returns>
        public T Determinant()
        {
            Debug.Assert(NonZeroCount == 0, "Zero Element Matrix");

            Number<T> Result = Number<T>.One;

            foreach (SparseElement se in this)
            {
                Debug.Assert(se.Column != se.Row, "Zero Element Matrix");

                Result *= (Number<T>)se.Value;
            }

            Debug.Assert(Number<T>.IsDefault(Result), "Matrix contains Zero's on the Diagonal!");

            return Result;
        }

        /// <summary>
        /// Returns the First Element of the Sparse Matrix.
        /// </summary>
        /// <param name="se">The First Element</param>
        /// <returns>False if there is no first element (NonZeroCount==0, after Clear()).</returns>
        public Boolean First(ref SparseElement se)
        {
            return First(ref se.Row, ref se.Column, ref se.Value);
        }

        /// <summary>
        /// Returns the First Element of the Sparse Matrix.
        /// </summary>
        /// <param name="r">The Row Index of the First Element</param>
        /// <param name="c">The Column Index of the First Element</param>
        /// <param name="v">The Value of the First Element</param>
        /// <returns>False if there is no first element (NonZeroCount==0, after Clear()).</returns>
        public Boolean First(ref Int32 r, ref Int32 c, ref T v)
        {
            this.FFN = 0;
            this.FFR = 0;

            if (FFN >= NonZeroCount)
            {
                return false;
            }

            v = Values[FFN];
            c = ColNdx[FFN];

            while (RowPtr[FFR] <= FFN)
            {                                                //= is sufficient!
                FFR++;
            }
            FFR--;

            r = FFR;

            return true;
        }

        /// <summary>
        /// Inverts the Matrix. 
        /// 
        /// Note: The Method used ONLY supports Diagonal Matrices without Zero values on the Diagonal.
        /// </summary>
        public void Invert()
        {
            foreach (SparseElement se in this)
            {
                Debug.Assert(se.Row == se.Column, "Matrix NOT Diagional!");

                //Determinant is not neccesary,
                //just replace A[ij] by 1/A[ij] for diagonal matrix
                //without zero's on the diagonal.

                Debug.Assert(Number<T>.IsDefault(se.Value), "Matrix contains Zero's on the Diagonal!");

                this[se.Row, se.Column] = Number<T>.Invert(se.Value);
            }
        }

        /// <summary>
        /// Returns the Next Element of the Sparse Matrix.
        /// </summary>
        /// <param name="se">The Next Element</param>
        /// <returns>False if there are no more elements.</returns>
        public Boolean Next(ref SparseElement se)
        {
            return Next(ref se.Row, ref se.Column, ref se.Value);
        }

        /// <summary>
        /// Returns the Next Element of the Sparse Matrix.
        /// </summary>
        /// <param name="r">The Row Index of the Next Element</param>
        /// <param name="c">The Column Index of the Next Element</param>
        /// <param name="v">The Value of the Next Element</param>
        /// <returns>False if there are no more elements.</returns>
        public Boolean Next(ref Int32 r, ref Int32 c, ref T v)
        {
            Boolean Result = false;

            FFN++;

            if (FFN < NonZeroCount)
            {
                v = Values[FFN];
                c = ColNdx[FFN];

                while (RowPtr[FFR] <= FFN)
                {
                    FFR++;
                }
                r = FFR - 1;

                Result = true;
            }
            return Result;
        }

        /// <summary>
        /// Removes Zero Values from a Sparse Array.
        /// 
        /// Note: We could do this in-place but that would be more complex.
        /// </summary>
        [Obsolete]
        public void CopyCompress()
        {
            Int32 cnt = 0;

            foreach (SparseElement se in this)
            {
                if (!Number<T>.IsDefault(se.Value))
                {
                    cnt++;
                }
            }

            //Only Compress when there are zero's!
            if (cnt != NonZeroCount)
            {
                //Use a copy of the Current Array....
                using (Sparse<T> copy = new Sparse<T>(this))
                {
                    this.Init(RowCount, ColCount, cnt, Comment);

                    foreach (SparseElement se in copy)
                    {
                        if (!Number<T>.IsDefault(se.Value))
                        {
                            this.AddElement(se);
                        }
                    }
                }
            }
        }

        private void Init(Int32 RowCount, Int32 ColCount, Int32 NonZeroCount)
        {
            Init(RowCount, ColCount, NonZeroCount, "");
        }

        private void Init(Int32 RowCount, Int32 ColCount, Int32 NonZeroCount, String Comment)
        {
            if (!Number<T>.IsNumeric)
            {
                throw new Exception("The type is not numeric");
            }

            this.ColCount = ColCount;
            this.RowCount = RowCount;
            this.NonZeroCount = NonZeroCount;
            this.Comment = Comment;

            this.Values = new T[NonZeroCount];
            this.ColNdx = new Int32[NonZeroCount];

            for (Int32 i = 0; i < NonZeroCount; i++)
            {
                Values[i] = default(T);
                ColNdx[i] = 0;
            }

            this.RowPtr = new Int32[RowCount + 1];

            for (Int32 i = 0; i < RowCount; i++)
            {
                RowPtr[i] = -1;
            }
            RowPtr[RowCount] = NonZeroCount;

            FNdx = 0;

            FFN = 0;
            FFR = 0;
        }

        /// <summary>
        /// Transposes the Matrix.
        /// </summary>
        public void Transpose()
        {
            //Use a copy of the Current Array....
            using (Sparse<T> copy = new Sparse<T>(this))
            {
                Init(copy.ColCount, copy.RowCount, copy.NonZeroCount, copy.Comment);

                for (Int32 r = 0; r < RowCount; r++)
                {
                    RowPtr[r] = 0;
                }

                //Count number nz in each row.
                for (Int32 r = 0; r < NonZeroCount; r++)
                {
                    RowPtr[copy.ColNdx[r]]++;                                           //zero based!
                }

                //Fill each cell with the starting point of the previous row.
                for (Int32 r = RowCount; r > 0; r--)
                {
                    RowPtr[r - 1] = RowPtr[r] - RowPtr[r - 1];
                }

                //Assign the new columns and values.
                Int32 ndx = 0;
                for (Int32 c = 0; c < copy.RowCount; c++)
                {
                    //for (Int32 ndx = copy.RowPtr[c]; ndx < copy.RowPtr[c + 1]; ndx++)
                    //TODO Look into this Dirty fix (ndx<nonZeroCount!)
                    while (ndx < NonZeroCount && ndx < copy.RowPtr[c + 1])
                    {
                        {
                            Int32 r = copy.ColNdx[ndx];
                            Int32 tr = RowPtr[r];

                            RowPtr[r]++;

                            ColNdx[tr] = c;
                            Values[tr] = copy.Values[ndx];

                            ndx++;
                        }
                    }
                }

                for (Int32 i = RowCount - 1; i > 0; i--)
                {
                    RowPtr[i] -= (RowPtr[i] - RowPtr[i - 1]);
                }
                RowPtr[0] = 0;
            }
        }

        /// <summary>
        /// Save Matrix as Harwell Boeing Sparse File.
        /// </summary>
        /// <param name="FileName">The FileName</param>
        public void SaveHarwellBoeing(String FileName)
        {
            using (new Profiling())
            {
                //Somehow we ended up with saving the other way round.
                Transpose();

                using (FileStream fs = new FileStream(FileName, FileMode.Create, FileAccess.Write, FileShare.None, 4 * 32768))
                {
                    StreamWriter sw = new StreamWriter(fs);

                    Int32 l1 = (Int32)Math.Ceiling(((1.0 * RowCount) + 1) / HBSplit);
                    Int32 l2 = (Int32)Math.Ceiling((1.0 * NonZeroCount) / HBSplit);
                    Int32 l3 = (Int32)Math.Ceiling((1.0 * NonZeroCount) / HBSplit);
                    Int32 l0 = l1 + l2 + l3;

                    //Header

                    //HB Matrix Type Identifier:
                    //
                    //First Character:
                    //      >R - Real matrix,
                    //       C - Complex matrix,
                    //       P - Pattern only (no numerical values supplied)
                    //Second Character:
                    //       S - Symmetric,
                    //       U - Unsymmetric,
                    //       H - Hermitian,
                    //       Z - Skew symmetric,
                    //      >R - Rectangular
                    //Third Character:
                    //      >A - Assembled
                    //       E - Elemental matrices(unassembled)

                    sw.WriteLine(("Harwell-Boeing: " + Comment).PadRight(72, ' ') + HBkey.PadRight(8));

                    sw.WriteLine(String.Format("{0} {1} {2} {3} {4}", l0.ToString().PadLeft(13), l1.ToString().PadLeft(13), l2.ToString().PadLeft(13), l3.ToString().PadLeft(13), 0.ToString().PadLeft(13)));
                    sw.WriteLine(String.Format("{0} {1} {2} {3} {4}", "rra".PadRight(13), ColCount.ToString().PadLeft(13), RowCount.ToString().PadLeft(13), NonZeroCount.ToString().PadLeft(13), 0.ToString().PadLeft(13)));
                    sw.WriteLine(String.Format("{0} {1} {2} {3}", "(8i)".PadLeft(15), "(8i)".PadLeft(15), "(8e)".PadLeft(19), "(8e)".PadLeft(19)));

                    //RowPtr
                    {
                        String s = "";
                        for (Int32 i = 0; i < RowCount + 1; i++)
                        {

                            s += RowPtr[i].ToString() + " ";
                            if ((i % HBSplit) == HBSplit - 1)
                            {
                                sw.WriteLine(s.Trim());

                                l1--;
                                s = "";
                            }
                        }

                        if (!String.IsNullOrEmpty(s))
                        {
                            sw.WriteLine(s.Trim());
                            l1--;
                        }
                    }

                    //ColNdx
                    {
                        String s = "";
                        for (Int32 i = 0; i < NonZeroCount; i++)
                        {
                            s += ColNdx[i].ToString() + " ";
                            if ((i % HBSplit) == HBSplit - 1)
                            {
                                sw.WriteLine(s.Trim());
                                l2--;
                                s = "";
                            }
                        }

                        if (!String.IsNullOrEmpty(s))
                        {
                            sw.WriteLine(s.Trim());
                            l2--;
                        }
                    }

                    //Values
                    {
                        String s = "";
                        for (Int32 i = 0; i < NonZeroCount; i++)
                        {
                            //s := s + Format(TSparseFmt + ' ', [FValues[i]]);
                            s += FormatSparseValue(Values[i]) + " ";
                            if ((i % HBSplit) == HBSplit - 1)
                            {
                                sw.WriteLine(s.Trim());
                                l3--;
                                s = "";
                            }
                        }

                        if (!String.IsNullOrEmpty(s))
                        {
                            sw.WriteLine(s.Trim());
                            l3--;
                        }
                    }

                    //Sanity Check
                    Debug.Assert((l1 == 0) && (l2 == 0) && (l3 == 0), "Harwell Boeing File Length Checks Failed.");

                    //AutoSave small matrices as NonSparse and CSV too. Always clean old ones on the fly...
                    String NonSparseFileName = Path.ChangeExtension(FileName, SparseExtension.NSext);
                    if (File.Exists(NonSparseFileName))
                    {
                        File.Delete(NonSparseFileName);
                    }

                    if ((RowCount * ColCount) < 1000)
                    {
                        SaveNonSparse(NonSparseFileName);
                        //SaveExcel(ChangeFileExt(outfile, CSext), rowlbls, collbls);
                    }

                    sw.Close();
                }

                //Somehow we ended up with saving the other way round, revert it.
                Transpose();
            }
        }

        private string FormatSparseValue(T value)
        {
            return String.Format(SparseFmt, value);
        }

        /// <summary>
        /// Save Matrix as MatrixMarket Sparse File.
        /// </summary>
        /// <param name="FileName">The FileName</param>
        public void SaveMatrixMarket(String FileName)
        {
            using (new Profiling())
            {
                using (FileStream fs = new FileStream(FileName, FileMode.Create, FileAccess.Write, FileShare.None, 4 * 32768))
                {
                    StreamWriter sw = new StreamWriter(fs);

                    sw.WriteLine(String.Format("%%{0}", Comment));
                    sw.WriteLine(String.Format("%%{0}", "G.W. van der Vegt 2010"));
                    sw.WriteLine(String.Format("{0} {1} {2}", RowCount, ColCount, NonZeroCount));

                    for (Int32 r = 0; r < RowCount - 1; r++)
                    {
                        for (Int32 j = RowPtr[r]; j < RowPtr[r + 1]; j++)
                        {
                            sw.WriteLine(String.Format("{0} {1} {2}", r, ColNdx[j], FormatSparseValue(Values[j])));
                        }
                    }

                    sw.Close();
                }
            }
        }

        /// <summary>
        /// Save the Matrix in a Non-Sparse Format.
        /// </summary>
        /// <param name="FileName">The FileName</param>
        public void SaveNonSparse(String FileName)
        {
            using (new Profiling())
            {
                //TODO Numeric Format.
                //TODO Correct for missing + sign (see decomposed matrices).
                //TODO Make tle line below match the SparseFmt width
                //TODO Make the SparseFmt width dynamic.
                //TODO Support Row/Column Lists.
                String nul = String.Format("      {0}", default(T));

                using (FileStream fs = new FileStream(FileName, FileMode.Create, FileAccess.Write, FileShare.None, 4 * 32768))
                {
                    StreamWriter sw = new StreamWriter(fs);

                    //TODO # of bits
                    sw.WriteLine(String.Format("'{0} (Sparseness={1}, Dimension={2}x{3} CxR, Floating Point Precision={4} bits)", Comment, SparseNess, ColCount, RowCount, 8));

                    String s = String.Format("{0:16}", @"[TERM\DOCUMENT]");
                    //if Assigned(files) and (files.Count = ColCount) then
                    //  for rr := 0 to Pred(ColCount) do
                    //    s := s + Format(' %*s', [TSparseFmtWidth + 4, ChangeFileExt(ExtractFileName(files.ValueFromIndex[rr]), '')])
                    //else
                    for (Int32 i = 0; i < ColCount; i++)
                    {
                        s += String.Format(" Col{0:000}", i);
                    }

                    sw.WriteLine(s);

                    Int32 r = 0;
                    Int32 c = 0;
                    T v = default(T);

                    First(ref r, ref c, ref v);
                    for (Int32 i = 0; i < RowCount; i++)
                    {
                        s = String.Format("Row{0:000}", i);

                        for (Int32 j = 0; j < ColCount; j++)
                        {
                            if ((i == r) && (j == c))
                            {
                                s += String.Format(" " + SparseFmt, v);
                                Next(ref r, ref c, ref v);
                            }
                            else
                            {
                                s += nul;
                            }
                        }

                        sw.WriteLine(s);
                    }

                    sw.Close();
                }
            }
        }

        /// <summary>
        /// Removes Zero Values from a Sparse Array.
        /// 
        /// Note: This method performs the Compression In-Place.
        /// </summary>
        public void Compress()
        {
            Int32 cvOffset = 0;
            Int32 rOffset = 0;

            //Three cases:
            //1) Nothing changes on a row.
            //2) We delete part of a row.
            //3) We delete the complete row.

            //Method:
            //1) We walk the RowPtr and enum the values of that row.
            //2) If an Element is zero we increment cvOffset and 
            //   copy all elements (and Colndx values) cvOffset to the left.
            //3) We count the number of values deleted and adjust upcoming 
            //   RowPtr's accordingly (rOffset).

            for (Int32 r = 0; r < RowCount; r++)
            {
                Int32 cnt = RowPtr[r + 1] - RowPtr[r];
                Int32 vals = RowPtr[r + 1] - RowPtr[r];

                for (Int32 ndx = RowPtr[r]; ndx < RowPtr[r + 1]; ndx++)
                {
                    if (Number<T>.IsDefault(Values[ndx]))
                    {
                        cvOffset++;
                    }
                    else
                    {
                        cnt--;

                        ColNdx[ndx - cvOffset] = ColNdx[ndx];
                        Values[ndx - cvOffset] = Values[ndx];
                    }
                }

                //Empty Row
                if (cnt == vals)
                {
                    //3) We deleted the complete Row
                }
                else if (cnt != 0)
                {
                    //2) We deleted some
                }
                else
                {
                    //1) Row Unchanged
                }

                //Adjust RowPtr
                //Note: The current RowPtr is not affected by the Deleted Values in it, 
                //      only the following rows. So we do this BEFORE adjusting rOffset so
                //      the incremented rOffset shows up during next Row.
                RowPtr[r] -= rOffset;

                rOffset += cnt;
            }

            //Adjust the NonzeroCount
            NonZeroCount -= cvOffset;

            //Adjust the internal Arrays.
            if (NonZeroCount == 0)
            {
                //We always leave at least one empty item in.
                NonZeroCount++;

                System.Array.Resize<Int32>(ref FColNdx, NonZeroCount);
                System.Array.Resize<T>(ref FValues, NonZeroCount);

                ColNdx[0] = 0;
                Values[0] = default(T);
            }
            else
            {
                System.Array.Resize<Int32>(ref FColNdx, NonZeroCount);
                System.Array.Resize<T>(ref FValues, NonZeroCount);
            }

            //Adjust RowPtr End marker!
            RowPtr[RowCount] = NonZeroCount;
        }

        /// <summary>
        /// Adds Zero Values to a Sparse Array in order to make it Non Sparse.
        /// </summary>
        public void Decompress()
        {
            using (Sparse<T> copy = new Sparse<T>(this))
            {
                //Create an Uncompressed Matrix.
                Init(copy.RowCount, copy.ColCount, copy.RowCount * copy.ColCount, copy.Comment);

                //Populate the non sparse Array so we can use Element[].
                for (Int32 r = 0; r < RowCount; r++)
                {
                    for (Int32 c = 0; c < ColCount; c++)
                    {
                        AddElement(r, c, default(T));
                    }
                }

                foreach (SparseElement se in copy)
                {
                    this[se.Row, se.Column] = se.Value;
                }
            }
        }

        /// <summary>
        /// CoFactors the matrix.
        /// 
        /// Note: The Method used ONLY supports Diagonal Matrices.
        /// </summary>
        public void Cofactor()
        {
            Debug.Assert(NonZeroCount == 0, "Zero Element Matrix");

            Number<T> dia = Number<T>.One;

            foreach (SparseElement se in this)
            {
                dia *= se.Value;
            }

            Debug.Assert(dia == (Number<T>)default(T), "Matrix contains Zero's on the Diagonal!");

            foreach (SparseElement se in this)
            {
                if (se.Column != se.Row)
                {
                    Debug.Assert(se.Column == se.Row, "Matrix NOT Diagional!");
                }

                //Divide by zero already caught by assertion on dia=0!
                this[se.Row, se.Column] = dia / se.Value;
            }

        }

        /// <summary>
        /// Fills this Matrix with the distances of a matrix. 
        /// </summary>
        /// <param name="a">Matrix a</param>
        /// <param name="aComment">Comment for the resulting Matrix</param>
        public void Distance(Sparse<T> a, String aComment)
        {
            //Use a copy of the Current Array....
            using (Sparse<T> copy = new Sparse<T>(a))
            {
                //Create a non-sparse Vector.
                Init(copy.RowCount, 1, copy.RowCount, aComment);

                //Populate the non sparse Array so we can use Element[].
                for (Int32 r = 0; r < RowCount; r++)
                {
                    AddElement(r, 0, default(T));
                }

                //Sum Sqr's 
                foreach (SparseElement se in copy)
                {
                    this[se.Row, 0] += ((Number<T>)se.Value * (Number<T>)se.Value);
                }

                //Take Sqrt's
                foreach (SparseElement se in this)
                {
                    Double sqrt = Math.Sqrt((Double)Convert.ChangeType(se.Value, typeof(Double)));

                    this[se.Row, 0] = (T)Convert.ChangeType(sqrt, typeof(T));
                }
            }
        }

        /// <summary>
        /// Fills this Matrix with the Dot Product of a matrix and a vector.
        /// </summary>
        /// <param name="a">Matrix a</param>
        /// <param name="b">Vector b</param>
        /// <param name="aComment">Comment for the resulting Matrix</param>
        public void Dot(Sparse<T> a, Sparse<T> b, String aComment)
        {
            Debug.Assert(b.RowCount == 1, "The Dot Product needs a Matrix and a Vector.");
            Debug.Assert(a.ColCount == b.ColCount, "The Dot Product needs a Matrix and a Vector of the same Length.");

            //Create a non-sparse Vector.
            Init(a.RowCount, 1, a.RowCount, aComment);

            //Populate the non sparse Array so we can use Element[].
            for (Int32 r = 0; r < RowCount; r++)
            {
                AddElement(r, 0, default(T));
            }

            //TODO Not VERY efficient...
            foreach (SparseElement sea in a)
            {
                foreach (SparseElement seb in b)
                {
                    if (sea.Column == seb.Column)
                    {
                        this[sea.Row, 0] += ((Number<T>)sea.Value * (Number<T>)seb.Value);
                    }
                }
            }
        }

        private Boolean InRange(Int32 ndx, Int32 min, Int32 max)
        {
            return (min <= ndx && ndx <= max);
        }

        /// <summary>
        /// Extracts a Rectangular area from an Array.
        /// 
        /// TODO: We could do an inplace version (set all except region to 0, compress and adjust internal arrays).
        /// </summary>
        /// <param name="a">The Array to Extract from</param>
        /// <param name="Left">Left</param>
        /// <param name="Top">Top</param>
        /// <param name="Right">Right</param>
        /// <param name="Bottom">Bottom</param>
        [Obsolete]
        public void Extract(Sparse<T> a, Int32 Left, Int32 Top, Int32 Right, Int32 Bottom)
        {
            //Use a copy of the Current Array....
            using (Sparse<T> copy = new Sparse<T>(a))
            {
                Right = Math.Min(Right, copy.ColCount);
                Bottom = Math.Min(Bottom, copy.RowCount);

                //Create a Non Sparse Matrix
                Int32 rows = Bottom - Top;
                Int32 cols = Right - Left;

                Init(rows, cols, rows * cols, copy.Comment);

                //Populate the non sparse Array so we can use Element[].
                for (Int32 r = 0; r < RowCount; r++)
                {
                    for (Int32 c = 0; c < ColCount; c++)
                    {
                        AddElement(r, c, default(T));
                    }
                }

                foreach (SparseElement se in copy)
                {
                    if (InRange(se.Row, Top, Bottom) &&
                        InRange(se.Column, Left, Right))
                    {
                        this[se.Row - Left, se.Column - Top] = se.Value;
                    }

                    if (se.Row > Bottom)
                    {
                        break;
                    }
                }

                //Make Matrix Sparse.
                Compress();
            }
        }

        /// <summary>
        /// Calculates the Frobenius Product bwteene a Matrix and a Vector.
        /// </summary>
        /// <param name="a">The Matrix</param>
        /// <param name="b">The Vector</param>
        /// <param name="aComment">Comment for the Resulting Matrix</param>
        public void Frobenius(Sparse<T> a, Sparse<T> b, String aComment)
        {
            Debug.Assert(b.ColCount == 1, "The Frobenius Product needs a Matrix and a Vector.");
            Debug.Assert(a.RowCount == b.RowCount, "The Frobenius Product needs a Matrix and a Vector of the same Length.");

            Distance(a, aComment);

            T fb = default(T);
            using (Sparse<T> tmp = new Sparse<T>())
            {
                tmp.Distance(b, "");

                Int32 r = 0;
                Int32 c = 0;
                tmp.First(ref r, ref c, ref fb);
            }

            foreach (SparseElement se in this)
            {
                this[se.Row, se.Column] = (Number<T>)se.Value / (Number<T>)fb;
            }
        }

        /// <summary>
        /// Fills this Matrix with the Cosine Similarity between a matrix and a vector.
        /// 
        /// NOTE: Untested
        /// </summary>
        /// <param name="a">Matrix a</param>
        /// <param name="b">Vector b</param>
        /// <param name="aComment">Comment for the resulting Matrix</param>
        public void Cosine(Sparse<T> a, Sparse<T> b, String aComment)
        {
            Debug.Assert(b.RowCount == 1, "The Cosine Similarity needs a Matrix and a Vector.");

            //TODO Look into this Assertion!
            //Debug.Assert(a.ColCount == b.ColCount, "The Cosine Similarity needs a Matrix and a Vector of the same Length.");

            //Vector
            Sparse<T> dp = new Sparse<T>();
            dp.Dot(a, b, "");

            //Vector
            Sparse<T> ul = new Sparse<T>();
            ul.Distance(a, "");

            //Number
            Int32 rr = 0;
            Int32 cc = 0;

            T qlv = default(T);
            Sparse<T> ql = new Sparse<T>();
            ql.Distance(b, "");
            ql.First(ref rr, ref cc, ref qlv);

            //Create a non-sparse Vector.
            Init(a.RowCount, 1, a.RowCount, aComment);

            //Populate the non sparse Array so we can use Element[].
            for (Int32 r = 0; r < RowCount; r++)
            {
                AddElement(r, 0, default(T));
            }

            foreach (SparseElement se in dp)
            {
                this[se.Row, 0] = (Number<T>)se.Value / ((Number<T>)ul[se.Row, 0] * (Number<T>)qlv);
            }
        }

        /// <summary>
        /// Multiplies two matrices.
        /// </summary>
        /// <param name="a">Matrix a</param>
        /// <param name="b">Matrix b</param>
        /// <param name="aComment">Comment for the resulting Matrix</param>
        public void Multiply(Sparse<T> a, Sparse<T> b, String aComment)
        {
            Debug.Assert(a.ColCount == b.RowCount, "The Arrays have incompatible Dimensions.");

            //Create an Uncompressed Matrix.
            Init(a.RowCount, b.ColCount, a.RowCount * b.ColCount, aComment);

            //Populate the non sparse Array so we can use Element[].
            for (Int32 r = 0; r < RowCount; r++)
            {
                for (Int32 c = 0; c < ColCount; c++)
                {
                    AddElement(r, c, default(T));
                }
            }

            b.Transpose();
            for (Int32 r = 0; r < RowCount; r++)
            {
                for (Int32 c = 0; c < ColCount; c++)
                {

                    Int32 k = a.RowPtr[r];
                    Int32 l = b.RowPtr[c];

                    while (l < b.RowPtr[c + 1] && k < a.RowPtr[r + 1])
                    {
                        if (a.ColNdx[k] == b.ColNdx[l])
                        {
                            Values[(r * ColCount) + c] += (Number<T>)a.Values[k] * (Number<T>)b.Values[l];

                            k++;
                            l++;
                        }
                        else if (a.ColNdx[k] < b.ColNdx[l])
                        {
                            k++;
                        }
                        else
                        {
                            l++;
                        }
                    }
                }
            }
            b.Transpose();

            Compress();
        }

        /// <summary>
        /// Extracts a Rectangular area from an Array.
        /// 
        /// Note: This method performs the Extraction In-Place.
        /// </summary>
        /// <param name="Left">Left</param>
        /// <param name="Top">Top</param>
        /// <param name="Right">Right</param>
        /// <param name="Bottom">Bottom</param>
        public void Reduce(Int32 Left, Int32 Top, Int32 Right, Int32 Bottom)
        {
            Int32 _Left = Left == -1 ? 0 : Math.Min(Left, ColCount);
            Int32 _Top = Top == -1 ? 0 : Math.Min(Top, RowCount);
            Int32 _Right = Right == -1 ? ColCount - 1 : Math.Min(Right, ColCount - 1);
            Int32 _Bottom = Top == -1 ? RowCount - 1 : Math.Min(Bottom, RowCount - 1);

            Console.WriteLine("Requested Reduce Edges: {0}x{1} .. {2}x{3}", Left, Top, Right, Bottom);
            Console.WriteLine("Actual Reduce Edges: {0}x{1} .. {2}x{3}", _Left, _Top, _Right, _Bottom);
            Console.WriteLine("Actual Reduce Dimension: {0}x{1}", _Right - _Left + 1, _Bottom - _Top + 1);

            Profiling.ConsoleStamp("Before Zeroing");

            //1 - Set all unneccesary values to zero.
            foreach (SparseElement se in this)
            {
                if (!(InRange(se.Row, _Top, _Bottom) && InRange(se.Column, _Left, _Right)))
                {
                    this[se.Row, se.Column] = default(T);
                }
            }

            Profiling.ConsoleStamp("After Zeroing, Before Compress");

            //2 - We need to compress to get the ColNdx and RowPtr arrays updated.
            Compress();

            Profiling.ConsoleStamp("After Compress, Before Leftmost");

            //3 - Clear Leftmost Columns
            if (_Left != 0)
            {
                for (Int32 i = 0; i < NonZeroCount; i++)
                {
                    ColNdx[i] -= _Left;
                }
            }

            Profiling.ConsoleStamp("After Leftmost, Before Rightmost");

            //4 - Clear Rightmost (Existing) Columns
            ColCount = _Right - _Left + 1;

            Profiling.ConsoleStamp("After Rightmost, Before Bottom");

            //5 - Clear Bottom Rows
            if (_Bottom <= RowCount - 1)
            {
                RowCount = _Bottom + 1;

                System.Array.Resize<Int32>(ref FRowPtr, RowCount + 1);
            }

            Profiling.ConsoleStamp("After Bottom, Before Top");

            //6 - Clear Top Rows
            while (RowCount > 0 && (RowPtr[0] == RowPtr[1]))
            {
                for (Int32 i = 1; i < RowPtr.Length; i++)
                {
                    RowPtr[i - 1] = RowPtr[i];
                }

                RowCount--;
                System.Array.Resize<Int32>(ref FRowPtr, RowCount + 1);
            }

            Profiling.ConsoleStamp("After Top");

            //TODO Sanity Check: NonZeroCount == RowPtr[RowCount];
            //TODO Sanity Check: Size ColNdx = NonZeroCount
            //TODO Sanity Check: Size Values = NonZeroCount
        }

        /// <summary>
        /// Save Matrix as CSV File.
        /// </summary>
        /// <param name="FileName">The FileName</param>
        public void SaveExcel(String FileName)
        {
            SaveExcel(FileName, null, null);
        }

        /// <summary>
        /// Save Matrix as CSV File.
        /// </summary>
        /// <param name="FileName">The FileName</param>
        /// <param name="rows">An optional list of Row Names</param>
        /// <param name="cols">An optional list of Column Names</param>
        public void SaveExcel(String FileName, StringCollection rows, StringCollection cols)
        {
            using (new Profiling())
            {
                using (FileStream fs = new FileStream(FileName, FileMode.Create, FileAccess.Write, FileShare.None, 4 * 32768))
                {
                    StreamWriter sw = new StreamWriter(fs);

                    String nul = String.Format(SparseFmt, default(T));

                    Int32 r = 0;
                    Int32 c = 0;
                    T v = default(T);

                    //Column Labels
                    String Result = @"[Row\Col]" + ExcelSeparator;
                    if (cols != null && cols.Count == RowCount)
                    {
                        foreach (String col in cols)
                        {
                            Result += EscapeCSV(col) + ExcelSeparator;
                        }
                    }
                    else
                    {
                        for (Int32 i = 0; i < ColCount; i++)
                        {
                            Result += String.Format("C{0}", i) + ExcelSeparator;
                        }
                    }

                    sw.WriteLine(Result.TrimEnd(new Char[] { ExcelSeparator }));

                    //Append Matrix
                    First(ref r, ref c, ref v);

                    for (Int32 i = 0; i < RowCount; i++)
                    {
                        Result = "";

                        //Row Labels
                        if (rows != null && rows.Count == RowCount)
                        {
                            Result += EscapeCSV(rows[i]);
                        }
                        else
                        {
                            Result += String.Format("R{0}", i);
                        }
                        Result += ExcelSeparator;

                        //Row Data
                        for (Int32 j = 0; j < ColCount; j++)
                        {
                            if ((i == r) && (j == c))
                            {
                                //We caught up with the NonZero Value.
                                Result += String.Format(SparseFmt, v);

                                Next(ref r, ref c, ref v);
                            }
                            else
                            {
                                //We still have a Zero Value.
                                Result += nul;
                            }

                            Result += ExcelSeparator;
                        }

                        sw.WriteLine(Result.TrimEnd(new Char[] { ExcelSeparator }));
                    }

                    sw.Close();
                }
            }
        }

        private String EscapeCSV(String value)
        {
            //See http://en.wikipedia.org/wiki/Comma-separated_values#Basic_Rules
            String Result = value;

            if (value.Contains('\"'))
            {
                Result = Result.Replace("\"", "\"\"");
            }

            if (value.Contains(ExcelSeparator) || value.Contains('"'))
            {
                Result = "\"" + Result + "\"";
            }

            return Result;
        }

        /// <summary>
        /// Dumps a Matrix as a String.
        /// </summary>
        /// <returns>The Matrix in String Format</returns>
        public override String ToString()
        {
            StringBuilder sb = new StringBuilder();

            Int32 r = 0;
            Int32 c = 0;
            T v = default(T);

            if (NonZeroCount != 0)
            {
                //Append Matrix
                First(ref r, ref c, ref v);

                for (Int32 i = 0; i < RowCount; i++)
                {
                    String Result = String.Format("{0:#######0}>", i);

                    for (Int32 j = 0; j < ColCount; j++)
                    {
                        if ((i == r) && (j == c))
                        {
                            Result += String.Format(" | {0:00.000}", v);
                            Next(ref r, ref c, ref v);
                        }
                        else
                        {
                            Result += String.Format(" | {0:00.000}", default(T));
                        }
                    }
                    Result += " |";

                    sb.AppendLine(Result);
                }

                //Append Statistics & Settings...
                sb.AppendLine(String.Format("{0} Columns", ColCount));
                sb.AppendLine(String.Format("{0} Rows", RowCount));
                if (NonZeroCount == 1 && Number<T>.IsDefault(Values[0]))
                {
                    sb.AppendLine(String.Format("{0} Zero Element Matrix", NonZeroCount));
                }
                else
                {
                    sb.AppendLine(String.Format("{0} NonZero Element Matrix", NonZeroCount));
                }
                sb.AppendLine(String.Format("Comment: '{0}'", Comment));
            }
            else
            {
                sb.AppendLine(String.Format("{0} Columns", ColCount));
                sb.AppendLine(String.Format("{0} Rows", RowCount));
                sb.AppendLine(String.Format("Zero Element Matrix"));
                sb.AppendLine(String.Format("Comment: '{0}'", Comment));
            }

            return sb.ToString();
        }

        /// <summary>
        /// Loads a rra Type HarwellBoeing Sparse Matrix.
        /// </summary>
        /// <param name="FileName">The filename</param>
        public void LoadHarwellBoeing(String FileName)
        {
            //0 ->Harwell-Boeing: DemoMatrix7                                             NetLsa  
            //1 ->            0             0             0             0             0 <<<WRONG, NO VALUES!
            //2 ->rra                       4             3             5             0
            //3 ->           (8i)            (8i)                (8e)                (8e)
            //4 ->0 1 3 5
            //5 ->3 1 3 1 3
            //6 ->5,0000 4,0000 -2,0000 -1,0000 11,0000

            String[] Lines = File.ReadAllLines(FileName);
            Char[] SplitOn = { ' ' };
            String[] Tokens = null;

            //Line 1
            Int32 l0 = 0;
            Int32 l1 = 0;
            Int32 l2 = 0;
            Int32 l3 = 0;

            //Line 2
            String type = "";
            Int32 Cols = 0;
            Int32 Rows = 0;
            Int32 Nzc = 0;

            //Line 3+
            Int32 ndx1 = 0;
            Int32 ndx2 = 0;
            Int32 ndx3 = 0;

            for (Int32 i = 0; i < Lines.Length; i++)
            {
                switch (i)
                {
                    case 0:                     //Header/Signature
                        //Skip
                        break;
                    case 1:                     //Lines.Length-3, l1, l2, l3, 0
                        Tokens = Lines[i].Split(SplitOn, StringSplitOptions.RemoveEmptyEntries);
                        if (Tokens.Length == 5)
                        {
                            l0 = Int32.Parse(Tokens[0]); //Lines.Length-3
                            l1 = Int32.Parse(Tokens[1]);
                            l2 = Int32.Parse(Tokens[2]);
                            l3 = Int32.Parse(Tokens[3]);
                        }
                        break;
                    case 2:                     //HBType, ColCount, RowCount, NonZeroCount, 0    
                        Tokens = Lines[i].Split(SplitOn, StringSplitOptions.RemoveEmptyEntries);

                        if (Tokens.Length == 5)
                        {
                            type = Tokens[0];   //HBType
                            Cols = Int32.Parse(Tokens[1]);
                            Rows = Int32.Parse(Tokens[2]);
                            Nzc = Int32.Parse(Tokens[3]);

                            Debug.Assert(type == "rra", String.Format("Unsupported Harwell Boeing Format: {0}.", type));

                            //Create Sparse Matrix.
                            Init(Rows, Cols, Nzc);
                        }
                        break;
                    case 3:                     //{8i..8e}
                        break;
                    default:                    //RowPtr, ColNdx, Values
                        Tokens = Lines[i].Split(SplitOn, StringSplitOptions.RemoveEmptyEntries);

                        //Sections never End/Start on one line.
                        if (ndx1 < RowCount + 1)
                        {
                            foreach (String Token in Tokens)
                            {
                                FRowPtr[ndx1++] = Int32.Parse(Token);
                            }
                        }
                        else if (ndx2 < NonZeroCount)
                        {
                            foreach (String Token in Tokens)
                            {
                                FColNdx[ndx2++] = Int32.Parse(Token);
                            }
                        }
                        else if (ndx3 < NonZeroCount)
                        {
                            foreach (String Token in Tokens)
                            {
                                FValues[ndx3++] = (T)Convert.ChangeType(Token, typeof(T));
                            }
                        }

                        break;
                }
            }

            Transpose();
        }

        /// <summary>
        /// Load an MatrixMarket Sparse Array.
        /// </summary>
        /// <param name="FileName">The filename to load</param>
        public void LoadMatrixMarket(String FileName)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Load an Excel compatible CSV File into a Sparse Array.
        /// 
        /// Note: This Method expects a Header Row and Column.
        /// </summary>
        /// <param name="FileName">The filename to load</param>
        public void LoadExcel(String FileName)
        {
            String[] Lines = File.ReadAllLines(FileName);

            //Extract RowCount from # of Lines, Compensatie for Header Row.
            Int32 RowCount = Lines.Length - 1;
            Int32 ColCount = 0;

            //Extract ColCount from Header, Compensate for Header Column.
            if (Lines.Length != 0)
            {
                ColCount = ParseCsvLine(Lines[0]).Length - 1;
            }
            else
            {
                //Console.WriteLine("Could not determine ColCount");

                return;
            }

            //Create a Non-Sparse Matrix.
            Init(RowCount, ColCount, RowCount * ColCount);

            //Process Rows containing Data, Skip Header Row.
            String[] Tokens = null;

            for (Int32 r = 1; r < Lines.Length; r++)
            {
                Tokens = ParseCsvLine(Lines[r]);

                //Process Columns, Skip Header Column.
                for (Int32 c = 1; c < Tokens.Length; c++)
                {
                    AddElement(r - 1, c - 1, (T)Convert.ChangeType(Tokens[c], typeof(T)));
                }
            }

            //Make Matrix Sparse.
            Compress();
        }

        /// <summary>
        /// Parses a single CSV line and tokenizes it correctly
        /// </summary>
        /// <param name="aLine">The CSV line to parse</param>
        /// <returns>The tokenized line</returns>
        private static String[] ParseCsvLine(String aLine)
        {
            Char[] Line = aLine.ToCharArray();

            List<String> Tokens = new List<String>();

            Int32 charpos = 0;
            Int32 cnt = 0;

            String Token = "";
            while (charpos < Line.Length)
            {
                switch (Line[charpos])
                {
                    case '\"':
                        {
                            Int32 quotes = 0;

                            while (charpos + 1 < Line.Length)
                            {
                                if (Line[charpos] == ';' && quotes % 2 == 0)
                                {
                                    //Add Token to List, Remove Enclosing Quotes First.
                                    Tokens.Add(Token.Substring(1, Token.Length - 2));

                                    Token = "";
                                    cnt++;

                                    charpos++;

                                    break;
                                }

                                //Process Quotes.
                                if (Line[charpos] == '\"')
                                {
                                    //And Escaped Quotes.
                                    if (Line[charpos + 1] == '\"')
                                    {
                                        //Skip 2nd \"
                                        charpos++;
                                    }
                                    else
                                    {
                                        quotes++;
                                    }
                                }

                                //Always add char
                                Token += Line[charpos];

                                //Advance charpos
                                charpos++;
                            }
                            break;
                        }

                    case ';':
                        {
                            //Add Token to List
                            Tokens.Add(Token);

                            Token = "";
                            cnt++;

                            charpos++;

                            break;
                        }

                    default:
                        Token += Line[charpos++];

                        break;
                }
            }

            if (!String.IsNullOrEmpty(Token))
            {
                Tokens.Add(Token);
                Token = "";
            }

            return Tokens.ToArray();
        }

        #endregion Methods

        #region AlgLib

        /// <summary>
        /// Load a 2 Dimensional Array into a Sparse Matrix.
        /// </summary>
        /// <param name="Value">The Array to load</param>
        public void SetAs2DArray(Double[,] Value)
        {
            Int32 nzcnt = 0;

            for (Int32 r = 0; r < Value.GetLength(0); r++)
            {
                for (Int32 c = 0; c < Value.GetLength(1); c++)
                {
                    if (!Number<T>.IsDefault((T)Convert.ChangeType(Value[r, c], typeof(T))))
                    {
                        nzcnt++;
                    }
                }
            }

            Init(Value.GetLength(0), Value.GetLength(1), nzcnt);

            for (Int32 r = 0; r < Value.GetLength(0); r++)
            {
                for (Int32 c = 0; c < Value.GetLength(1); c++)
                {
                    T value = (T)Convert.ChangeType(Value[r, c], typeof(T));
                    if (!Number<T>.IsDefault(value))
                    {
                        AddElement(r, c, value);
                    }
                }
            }
        }

        /// <summary>
        /// Load a 1 Dimensional Array into a Sparse Matrix.
        /// 
        /// NOTE: The resulting array will be a diagonal one.
        /// </summary>
        /// <param name="Value">The Array to load</param>
        public void SetAs1DArray(Double[] Value)
        {
            Int32 nzcnt = 0;

            for (Int32 i = 0; i < Value.Length; i++)
            {
                if (!Number<T>.IsDefault((T)Convert.ChangeType(Value[i], typeof(T))))
                {
                    nzcnt++;
                }
            }

            Init(Value.Length, Value.Length, nzcnt);

            for (Int32 i = 0; i < Value.Length; i++)
            {
                T value = (T)Convert.ChangeType(Value[i], typeof(T));
                if (!Number<T>.IsDefault(value))
                {
                    AddElement(i, i, value);
                }
            }
        }

        /// <summary>
        /// Return the Matrix as a 2 Dimensional Array.
        /// </summary>
        /// <returns>The 2 Dimensional Array</returns>
        public Double[,] GetAs2DArray()
        {
            Double[,] Result = new Double[RowCount, ColCount];

            for (Int32 r = 0; r < RowCount; r++)
            {
                for (Int32 c = 0; c < ColCount; c++)
                {
                    Result[r, c] = default(Double);
                }
            }

            SparseElement se = new SparseElement();

            if (First(ref se))
            {
                Result[se.Row, se.Column] = (Double)Convert.ChangeType(se.Value, typeof(Double));

                while (Next(ref se))
                {
                    Result[se.Row, se.Column] = (Double)Convert.ChangeType(se.Value, typeof(Double));
                }
            }

            //foreach (SparseElement se in this)
            //{
            //    Result[se.Row, se.Column] = (Double)Convert.ChangeType(se.Value, typeof(Double));
            //}

            return Result;
        }

        #endregion Alglib

        #region IEnumerator Members

        /// <summary>
        /// Returns an Enumerator for this Object. Due to a quirk in the implementation of IEnumerate 
        /// we have to call Reset() manually before returning the Enumerator.
        /// </summary>
        /// <returns>The IEnumerator of this object</returns>
        public IEnumerator GetEnumerator()
        {
            try
            {
                //If implemented in the same class there can be only one instance. 
                //So we have to call Reset manually or a 2nd foreach will not return any elements.
                Reset();

                return (IEnumerator)this;
            }
            finally
            {
                // cleanup code here 
            }
        }

        /// <summary>
        /// Contains a Sparse Element from the Sparse Matrix.
        /// </summary>
        public struct SparseElement
        {
            internal Boolean Reset;

            /// <summary>
            /// The Row Index of the Sparse Element
            /// </summary>
            public Int32 Row;

            /// <summary>
            /// The Column Index of the Sparse Element
            /// </summary>
            public Int32 Column;

            /// <summary>
            /// The Value of the Sparse Element
            /// </summary>
            public T Value;
        }

        private SparseElement fCurrent = new SparseElement();

        /// <summary>
        /// Returns the Current SparseElement
        /// </summary>
        public object Current
        {
            get { return fCurrent; }
        }

        /// <summary>
        /// Moves to the next SparseElement
        /// </summary>
        /// <returns>True if there was a next SparseElement</returns>
        public bool MoveNext()
        {
            if (fCurrent.Reset == default(Boolean))
            {
                fCurrent.Reset = !fCurrent.Reset;
                return First(ref fCurrent.Row, ref fCurrent.Column, ref fCurrent.Value);
            }
            else
            {
                return Next(ref fCurrent.Row, ref fCurrent.Column, ref fCurrent.Value);
            }
        }

        /// <summary>
        /// Resets the Enumerators
        /// </summary>
        public void Reset()
        {
            fCurrent.Reset = default(Boolean);
        }

        /// <summary>
        /// Just a Dummy Dispose for now.
        /// </summary>
        public void Dispose()
        {
            //TODO
        }

        #endregion

        #region Other

        //There are 3 Parse Methods, we need Parse(String)!
        //MethodInfo miT = typeof(T).GetMethod("Parse",
        //                                    BindingFlags.InvokeMethod | BindingFlags.Public | BindingFlags.Static,
        //                                    null,
        //                                    CallingConventions.Any,
        //                                    new Type[] { typeof(String) },
        //                                    null);
        //Object[] parms = { sqrt.ToString() };
        //T NewValue = (T)miT.Invoke(null, parms);

        #endregion Other

        #region IDisposable Members

        void IDisposable.Dispose()
        {
            this.Dispose();
        }

        #endregion
    }

    public static class SparseExtension
    {
        public static readonly String HBext = ".hb";
        public static readonly String NSext = ".ns";
        public static readonly String MMext = ".mm";
        public static readonly String CSVext = ".cvs";
        public static readonly String SLext = ".ls";
        public static readonly String LSext = ".lsa";
    }

}
